      !< @brief Rotina de suavização do modelo de velocidades
      !> O valor da velocidade em um nó é calculado como a média da velocidade dos vizinhos.
      !> Este processo é reptido nsmooth vezes.
      subroutine suaviza(c,nx,nz)
      implicit none

      double precision:: c(nx,nz)
      double precision, dimension(:,:), allocatable:: caux
      integer :: nx,nz
      integer :: i,j,ksmooth
      integer :: ierr !< Indica se houve erro na alocação de memória
      integer, parameter ::  Nsmooth=2
 
      allocate(caux(nx,nz),stat=ierr)


      do ksmooth = 1, Nsmooth
	  j=1
	      do i=1,nx
		  if(i.eq.1)then
		      caux(i,j) =( c(i+1,j)  + c(i,j+1))/2.0d0
		  elseif (i.eq.nx)then
		      caux(i,j) =( c(i-1,j)  + c(i,j+1))/2.0d0
		  else
		      caux(i,j) =( c(i-1,j) + c(i+1,j) + c(i,j+1))/3.0d0
		  endif
	      enddo

	  do j=2,nz-1
	      do i=1,nx
		  if(i.eq.1)then
		      caux(i,j) =( c(i+1,j) + c(i,j-1) + c(i,j+1))/3.0d0
		  elseif (i.eq.nx)then
		      caux(i,j) =( c(i-1,j) + c(i,j-1) + c(i,j+1))/3.0d0
		  else
		      caux(i,j) =( c(i-1,j) + c(i+1,j) + c(i,j-1) + c(i,j+1))/4.0d0
		  endif
	      enddo
	  enddo

	  j=nz
	      do i=1,nx
		  if(i.eq.1)then
		      caux(i,j) =( c(i+1,j)  + c(i,j-1))/2.0d0
		  elseif (i.eq.nx)then
		      caux(i,j) =( c(i-1,j)  + c(i,j-1))/2.0d0
		  else
		      caux(i,j) =( c(i-1,j) + c(i+1,j) + c(i,j-1))/3.0d0
		  endif
	      enddo

	  do j=1,nz
	      do i=1,nx
		  c(i,j) = caux(i,j) !< atualiza o modelo
	      enddo
	  enddo

      enddo


      deallocate(caux)

      return

      end subroutine suaviza

